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Stochastic fluctuations are central to the understanding of extinction dynamics. In the context of 
population models they allow for the description of the transition from the vicinity of a non-trivial 
fixed point of the deterministic dynamics to a trivial fixed point, where the population has become 
extinct. To characterize analytically the fluctuations of a given stochastic population model, one 
can operate within the so-called linear-noise approximation. Here the fluctuations are taken to be 
Gaussian, and the phenomenon of extinction is so rare as to be negligible, for all practical purposes. 

When the size of the population becomes small, non-Gaussian fluctuations are instead found. Two 
analytical schemes are in principle available to determine the nature of the distribution of associated 
fluctuations: the generalized van Kampen system-size expansion beyond the conventional order of 
approximation and a WKB-like method. Here we investigate the accuracy of these two different 
approximation schemes, with reference to a simple stochastic process that converges to the logistic 
equation in the deterministic limit. 
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I. INTRODUCTION 

The huge increase recently seen in the construc¬ 
tion of quantitative models in the biological sciences, 
especially ecology, has prompted a renewed interest 
in stochastic effects in these systems flj. In com¬ 
puter simulations processes such as birth, death, pre¬ 
dation, and many others, are quite naturally treated 
as random. In the mathematical versions of these 
individual-based models, the stochastic effects are 
present; they only vanish in the limit of infinitely 
many individuals, when deterministic equations are 
obtained. However, in some important situations, the 
stochastic effects can be so large as to invalidate the 
conclusions based on a purely deterministic analy¬ 
sis SH3- Hence the interest in stochastic dynamics, 
and the mathematical techniques which can be used 
to elucidate them. 

The deterministic limit usually yields a set of cou¬ 
pled ordinary differential equations in time. These 
can be analysed using the well-developed ideas from 
the theory of dynamical systems involving concepts 
such as trajectories in phase space and attractors of 
various kinds. Stochastic systems have their own, 
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different, set of concepts such as stochastic fluctua¬ 
tions about the deterministic trajectories, rare events 
which lead to transitions from one attractor of the 
deterministic dynamics to another, stationary prob¬ 
ability distributions, and so on. It is the tools that 
can be used to calculate some of these quantities for 
a wide class of models that will of interest to us in 
this paper. 

We will focus on extinction dynamics, that is, tran¬ 
sitions from the vicinity of a non-trivial fixed point 
of the deterministic dynamics, to a trivial fixed point 
where all the individuals have become extinct. If the 
parameters of the models are chosen so that the two 
fixed points are not too close, then extinctions will 
be rare and the state of the system will continue to 
fluctuate about the non-trivial fixed point for a long 
period of time. In this situation the van Kampen 
system-size expansion Q] has proved to be a powerful 
tool. Essentially the method gives the deterministic 
equations to leading order in a expansion in (inverse) 
system-size, and linear stochastic corrections to this 
result at next-to-leading order. This linear approxi¬ 
mation, sometimes called the linear-noise approxima¬ 
tion, corresponds to Gaussian fluctuations, and is an 
excellent approximation if extinctions are so rare as 
to be negligible. 

In principle, extinction events can be incorporated 
into the system-size expansion by going beyond the 
Gaussian approximation, and including non-linear 
terms in the stochastic differential equation giving 
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the approximate stochastic dynamics. This should 
give rise to a tail on one side of the probability dis¬ 
tribution function which will characterise the extinc¬ 
tion process. These higher-order calculations have 
only been carried out recently, and then only by a 
few authors &S3- In addition if there are only a 
few individuals in the system, extinction effects will 
be very important, and one may have to go to quite 
high-order in the expansion to get an accurate form 
for what will be a very non-Gaussian distribution. 

By contrast, the standard technique to look at rare 
events, found going from one metastable state to an¬ 
other, is to use a WKB-like approximation. This 
goes under many names: large deviation theory (l4l j , 
the instanton method (l5j, Freidlin-Wentzell theory 
0, amongst others. It consists of postulating that 
the dominant contribution to the probability distri¬ 
bution is exponentially small in IV, the system-size, 
that is, is of the order of e~ NS . Here S is a function 
of variables describing the system in the determinis¬ 
tic limit which turns out to satisfy a Hamilton-Jacobi 
equation. The corresponding Hamiltonian, known as 
the FreidlinWentzell Hamiltonian, can be used to de¬ 
scribe the extinction trajectories, even though they 
are stochastic in nature. 

The aim of this paper is to explore the connec¬ 
tion between the van Kampen system-size expansion 
at higher-order and the WKB method. They have 
a very different basis, and to the best of our knowl¬ 
edge, their predictive power in regimes where extinc¬ 
tions are important have not been compared. We will 
carry out the explicit assessment of their range of va¬ 
lidity and comparison with numerical simulations on 
a specific stochastic system with one degree of free¬ 
dom to minimise numerical errors. In the case of the 
WKB method, most of the steps can be performed 
analytically, which also makes the interpretation of 
the results more straightforward. A simple dynami¬ 
cal system with a stable non-trivial fixed point and 
an unstable trivial fixed point is the logistic equation 
<j> = rcf> [1 — (<j>/k)\, and the stochastic model we will 
choose to study will have this equation as its deter¬ 
ministic limit. 


II. MODEL 

In this section we will introduce the individual 
based model we will analyse and write down the mas¬ 
ter equation, which governs its stochastic dynamics. 
The two techniques we are comparing can be viewed 
as different approximations to these dynamics. We 
will describe them in turn and then compare them to 
numerical simulations of the original individual based 
model in Section EU 

Following the discussion in the Introduction, prob¬ 
ably the simplest model which contains the features 
which we wish to explore is a system containing iden¬ 


tical individuals, which we will denote by A. We as¬ 
sume that there are n such individuals, and since the 
size of the system is taken to be characterised by an 
integer N, we suppose that there are (N — n) nulls 
denoted by E. These are vacancies, which in a spatial 
version of the model would denote spaces which could 
potentially be colonised by a individual. If the only 
processes are (asexual) birth, competition and death, 
then we may define the well-mixed model through the 
reactions 

A + E A A + A, A + A A- A + E, aAe. ( 1 ) 

The last equation, for instance, indicates that an in¬ 
dividual of type A dies at a rate d to give a vacancy, 
E. Simple combinatorics then gives the rate at which 
the number of individuals increases from n to (n+ 1) 
to be given by b ( n/N ) (N — n)/N. A more accurate 
statement of these rates would replace one of the N 
factors in the denominator by (N — 1), but since we 
wish to keep the analysis as simple as possible we will 
not do this. If we scale the time by a factor of N, then 
the transition rate from state n to state (n + 1) may 
be written as: 

T{n+l\n) = bn(l- (2) 

Similarly, the transition rate from state n to state 
(n — 1) is 


T(n- l\n) = n (d + c^j , ( 3 ) 


where once again factors of N(N — 1) and n(n — 1) 
have been replaced by N 2 and rr respectively. Since 
the transition rates T(n+l\n) and T(n—l\n) define 
the model, this choice simply corresponds to a slight 
variant of the standard model, which can be justified 
on grounds of simplicity. 

The master equation is an equation for the rate of 
change with time of the probability of finding n indi¬ 
viduals in the system at time t, denoted by P(n,t)- 
Since this is simply the rate of transitions into the 
state n minus the rate of transitions out of state n, it 
reads @ 


dP(n, t) 
dt 


T(n\n + 1 )P(n + 1, t) 

T{n\n — l)P(n — 1, t) (4) 

T(n — 1| n)P(n, t ) — T(n + 1| n)P(n, t). 


This equation cannot be solved exactly, so we need to 
resort to either numerical methods or approximation 
techniques. It is frequently simpler to simulate 0] 
the processes given in Eq. ©, rather than numeri¬ 
cally solve the master equation, and the results we 
give in Section [TTH to assess the accuracy of the ap¬ 
proximation techniques are found in this way. We 
now briefly outline the two approximation methods 
that we use in this paper. 
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A. The van Kampen system-size expansion 


The first, the van Kampen system-size expansion, 
has as the leading-order approximation the determin¬ 
istic differential equation found by taking the limit 
N —> oo. In the case of the model just described, 
this is the logistic equation given in the Introduc¬ 
tion. However this equation emerges as the leading 
order approximation to the model defined by Eqs. © 
and ©, and does not have to be postulated inde¬ 
pendently. The next-to-leading order gives a lin¬ 
ear stochastic differential equation, which describes 
Gaussian fluctuations about the deterministic result. 
If the intention is to simply study the stochastic dy¬ 
namics of the model well away from the boundaries, 
it is usually sufficient to work to this order. How¬ 
ever it is possible to go higher orders to obtain non- 
Gaussian corrections to the probability distribution 
function (pdf). One of the main aims of this paper 
is to argue that these higher-order corrections enable 
reliable estimates for the pdf to be obtained very close 
to the boundaries. 

To apply the van Kampen expansion we first write 
down the master equation © in terms of step- 
operators £ ± defined by £ ± f(n ) = /(n± 1 ), where / 
is an arbitrary function Q, 

dP ^ t ’ ^ = ( £+ ~ X ) t T ( n _ X l n) p (n, *)] 

+ (£~ — l) [T(n + l\n)P(n,t)]. (5) 


The ansatz which forms the basis of the method is to 
write Q 


— = 0(i) + -^L. 


( 6 ) 


Here 0(f) is the solution of the deterministic equation 
valid in the limit TV — > oo and £ is the (continuous) 
stochastic variable which gives the deviation of the 
stochastic trajectory from this deterministic value. 
The pdf when written in terms of £ is denoted as 
n(£, £), thus P(n,t) = n(£, t). After the change of 
variables the left-hand side of the master equa¬ 
tion © becomes: 


dP an r— an d0 

dt dt a£ dt' 


( 7 ) 


The right-hand side of the master equation in the 
form © can also be written in terms of 0 and £ by 
(i) eliminating n in the transition rates © and © 
using Eq. ©. and (ii) noting that the step-operators 
may be written as 


£± 


^ (±i)* i a f 
+ Z^ £! N ( / 2 dP' 


( 8 ) 


Equating the left-hand and right-hand sides of the 
master equation, after rescaling time by introducing 


t = t/N, we can match inverse powers of TV 1 / 2 to 
obtain a set of equations for the dynamics of the pro¬ 
cess. This is carried out explicitly in Appendix [A] At 
leading order — obtained by matching the coefficients 
of TV -1 / 2 — one finds the equation: 


d(j) 

dr 



(9) 


where r = b — d and K = (6 — d)/{b + c). This is 
the logistic equation, which could have been guessed 
as the deterministic limit of the model, even if the 
identification of the constant K is not so obvious. 
This has the solution: 




[I\ - 0 q] e~ rT + 0 O ’ 


( 10 ) 


where 0o = 0(0). This has the required feature that, 
as long as 0 o 7 ^ 0 , then 0 (r) —> 0 * as t —> 00 , where 
0* = K is the non-trivial fixed point. 

Once the leading-order contributions have been ex¬ 
tracted the left-hand side is simply dH/dr, but the 
right-hand side contains derivatives of n with respect 
to £ of all orders. The resulting equation has the gen¬ 
eral structure: 


5n V 1 1 <9 fe+1 n 

dr ~ ^ (k+ 1 )! JV(fc-i)/ 2 -' fc+1 ^ d £k +1 

00 ' 1 1 gk 

^ ^2 T\ N(k- 1)/2^(^) QCk 
k= 1 ’ S 

00 1 1 afc—1 

+ ^ (fc - 1)! g^k -1 n ] >( n ) 


where the explicit form of the functions /&, gt and the 
constant hk are given in Appendix [A] 

To proceed in the analysis we introduce the mo¬ 
ment of order q of the distribution n: 


<£*> = yVn(£)d£. (12) 


From the generalized Fokker-Planck equation ED 
one can obtain a set of ordinary differential equations 
for the coupled evolution of the moments of the dis¬ 
tribution n. The method is straightforward @] and 
consists in multiplying both sides of Eq. ED by £« 
and integrating by parts over the variable £. One 
finds 


d(€ q ) (—l) fc+ 1 g !A. + 1 (0)(£^ fc+1 >) 

dr E (fc + l)!AT(fc- 1 )/2( g _(fc + i))! 
(-1 ) fc g!fffc(£ 9 ~ fc+1 ) 

E k \ N (k--L)/2(q- k y 

^ (-1 ) fc - 1 g !/ lfc _ 1 (£«-( fe - 1 )+ 2 ) 

+ ^ (k- 1 )\(q -{k- 1 ))!iV( fc - 1 )/ 2 ' 


( 13 ) 
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It is clear from Eq. m, that the equation for 
d(^ q )/d,T depends on (£ 9+1 ) and so the system of 
equations does not close at any finite value of q. The 
general way to proceed in such cases is to impose 
some form of truncation which will eventually lead to 
a closed, self-consistent set of equations. For the case 
of the model we are investigating here, the fact that 
Eq. (fT51) is an expansion in powers of 1/y/N, sug¬ 
gests the possibility of implementing a rather natural 
truncation scheme, which we will return to in Sec¬ 
tion m Once complemented by a particular closure 
scheme, the system of Eqs. USD, up to a given value 
of q , can be solved numerically and the estimated mo¬ 
ments used to reconstruct the distribution via Fourier 
inversion, as discussed in Section m 


B. The WKB expansion 

The second method is to use a WKB approxima¬ 
tion. This approximation is well-known in the study 
of ordinary differential equations [l8j, but has also 
been extensively used in the analysis of rare events in 
stochastic systems, both by mathematicians 16], and 
theoretical physicists using path-integrals Pi|. Here 
will look at a variant of the method which starts off 
from the discrete state master equation, rather than 
from the Fokker-Planck equation and will follow the 
methodology described in Refs. [lf| and |20J. 

The approximation involves a different form of scal¬ 
ing to that described in Section Hi Al Now, the start¬ 
ing point is the master equation ([5]), but written in 
terms of x = n/N and N, rather than n. Specif¬ 
ically, we write P(n,t) = P(Nx,t) = n(x,i) and 
T±(n ± l|n) = NQ±(x), so that: 

f2+(x) = bx(l — x), D_(x) = x(d + cx). (14) 


respect to N 1 and collecting together the leading 
order terms yields 

'Y, f \(x) [exp(r5' , (x) — 1)] = 0, (17) 

r=z Ll 


where S'(x) = dS(x)/dx. The above equation can 
be seen as a stationary Hamilton-Jacobi equation 
H{x, S"(x)) = 0 for an action S with Hamiltonian 

H(x,p)= Y ftr(z)[exp(rp) - 1], (18) 

r=±l 


where p = S' (x). FromEq. m one gets the following 
Hamilton’s equations: 


x 

P 


dH _ 
dp 

dH 

dx 


Y^ r fl r (x) exp(rp), 

r=± 1 

= - Y [ ex P ( r P) ~ !] 

r=±l 


dtl r (x) 

dx 


, ( 19 ) 


where the dot denotes differentiation with respect to 
time. 

From the solution of these equations one finds the 
so-called fluctuation trajectories Xf and the corre¬ 
sponding momenta p /. For the zero-energy solution 
H = 0 that we are interested in, the action calcu¬ 
lated along a given fluctuation trajectory, starting at 
the fixed point at time to, is given by 

S f = [ pfXfdt', (20) 

Jto 

and so can be calculated from a knowledge of pf and 
*/• 

From Hamilton’s equations uni) one can readily 
check that a trivial solution, Pf = 0 exists, provided 
that 


We assume that N is sufficiently large that x is effec¬ 
tively continuous. Looking for the quasi-stationary 
solution of the master equation m, in the basin of 
the attraction of the stable fixed point, leads to the 
following equation: 

a - (* + v) 1n (* + v) + n+ ( x ~ v) 1n ( x - v) 

— [n_(x) + n+(x)] n(x) = o. 

(15) 


To solve the above time-independent equation for 
the stationary distribution, we apply the WKB ap¬ 
proximation 119), 20] by assuming the following form 
for n(x): 

n(x) = K(x)exp(-NS(x)) , (16) 


where both S(x) and K (x) are of the order of unity. 
Substituting (USD into m, Taylor expanding with 


x = f2+(x) — f2_(x), (21) 

where from now on we drop the subscript /. 

This is customarily called the relaxation trajectory 
and corresponds to the deterministic (N —> oo) ap¬ 
proximation of the model. In fact, multiplying both 
sides of the original master equation by n and sum¬ 
ming over all possible states one finds that 

^ = JVn+ (x) -Nil- (x), (22) 

dt 

where (n) = i)- Dividing then by N and 

taking the limit N —> oo, one finds Eq. (1 21 [) . This 
equation is nothing else but the logistic equation, and 
the relaxation trajectory eventually converges to the 
stable fixed point. 

In this paper we are interested rather in the so¬ 
lution of the Hamilton’s equations (HU) with p / 0. 
This enables us to explore trajectories which not al¬ 
lowed in the deterministic limit. This solution, which 








5 


has a non-trivial value of the momentum, can be 
found by setting H = 0 in Eq. m which leads to a 
quadratic equation in e p . Solving this equation yields 

P = In o ; z = ^-(®) - 0+(®)- (23) 

U+[x) 

The actions of these paths can now be calculated. 
In the case p = 0 it is simply zero, and in the case 


given by Eq. (TZ!Hl 

one finds that 


S(x) — S(x*) = 

(d + cx) ln(d + cx) 

(d + cx) 

c 

c 

+ 

(b — bx ) ln(6 — bx) 

(b — bx) 

b 

b 

+ 

{d + c) (c + d) 

h 111 

c b 

b(d + c) 
c+b 


(24) 

Here we have used the expression x* = (b — d)/(b + c) 
for the non-trivial fixed point, found from the condi¬ 
tion S2 + (x*) = fl_(x*). 

The next-to-leading terms in the WKB approxima¬ 
tion (flbl) is the prefactor K(x). It is found pH [20] to 
obey the equation 

dH K' _ 1 ,d 2 H d 2 H 
dp K dp 2 dpdx 


and therefore: 

n (z) ~ exp[-iVS(a:*)] exp ( - ^(x - x*) 2 ), 

(29) 

where A = N(b + c) 2 /b(c + d). By imposing the nor¬ 
malization condition for n(x) we obtain the following 
expression for A: 


A = \j — exp [NS(x*)] fi+(x*), 


which eventually yields 


(30) 


1NS"{x*) 

[H + (x*)H_(x*)i 

V 2t r 

. H+(x)fl_(x) . 


n(x) = 


x exp(— N[S(x) — S(x*)]). 


1/2 


(31) 


In conclusion, we have derived a closed analytical 
expression for the quasi-stationary distribution n, us¬ 
ing the WKB approximation procedure. In the next 
section we will test the adequacy of formula m , as 
well as the corresponding result obtained within the 
framework of the van Kampen system-size expansion, 
by performing a direct comparison with the outcome 
of stochastic simulations. 


III. NUMERICAL SIMULATIONS 


For the case of interest here, p ^ 0, with p(x) speci¬ 
fied by equation (1251) . one can solve equation (05l) to 
obtain: 


K (x) 


A 

( x ) ’ 


(26) 


where A is a constant to be determined. For the sake 
of completeness we also give the expression for K(x) 
which applies to the case p = 0. It is 


K(x) 


B 

f2+(x) — fl-(x) 1 


(27) 


where B is another constant. 

We can find the constant A by normalising the 
probability distribution. If we expand the quasi- 
stationary distribution n(x) on the trajectory with 
p / 0 and in the vicinity of the non-trivial fixed point 
x = x* we have that 


n (z) 


N 


^/Q+(x*)^LAx*) 


exp 


- NS(x*) 


(28) 


- -S"{x*)(x-x*) 2 


where use had been made of the fact that S'{x*) = 
ln(fi_(x*)/fi+(x*)) = 0. A straightforward calcula¬ 
tion allows us to write: 


S"(x*) = (b + c) 2 /b{c + d ), 


In Section |H] by extending the van Kampen 
system-size expansion beyond the Gaussian order, we 
obtained a set of coupled differential equations for 
the moments of the distribution of fluctuations. The 
knowledge of these moments enables us in principle 
to reconstruct the corresponding distribution via a 
standard Fourier inversion. In this Section we will 
compare this theoretical prediction, along with the 
WKB result (13TT) . and assess their validity through 
direct simulation of the stochastic dynamics. 

While it is straightforward to display the result of 
the WKB analysis, some comments are in order con¬ 
cerning the interpretation of the calculation based on 
the van Kampen expansion. To find the stationary 
distribution of fluctuations we consider the ensemble 
of the first q moments, and impose a truncation in the 
van Kampen expansion by omitting the term propor¬ 
tional to (£ q+1 )/VN in the equation for d(^ q )/dT. 
This is the k = 2 contribution in the last sum on the 
right-hand side of equation ra¬ 
in principle, one cannot formally drop such a term, 
but we can assess its importance since we are work¬ 
ing within a large system size (large N ) approxima¬ 
tion. It is found that the structure of the governing 
equations means that the effect of neglecting a con¬ 
tribution proportional to 1 /y/N in the equation for 
the q -th moment, becomes rapidly less important as 
the order of the moments decreases. So by taking q 
large enough, we would expect that the errors made 
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in the estimates of a significant number of the lowest 
order moments would be negligible. However, these 
arguments are heuristic, and ultimately only a com¬ 
parison with simulations will provide an a posteriori 
validation of the proposed approximation. 

As an additional remark, we recall that we are 
interested here in the the asymptotic distribution 
of fluctuations around the non-trivial fixed point of 
the deterministic dynamics. This implies setting the 
derivatives d{^ q )/dr in Eq. (fl3l) to zero, after hav¬ 
ing imposed the fixed point condition <j> = cf>* on 
the functions fk and g defined in Appendix A. 
The system of differential equations for the evolution 
of the moments is hence transformed into an alge¬ 
braic system that can be readily solved by matrix 
inversion. The distribution of fluctuations is finally 
determined by Fourier inverting the corresponding 
moment-expansion. In the following, we will report 
results obtained when considering the first thirty mo¬ 
ments in the expansion, i.e., taking the first q = 30 
algebraic equations for the stationary moments. 

In Fig. [Tj we compare the distributions of fluctu¬ 
ations obtained via the van Kampen procedure, at 
different order of approximations of Eq. (1131) . The 
usual van Kampen system-size expansion assumes a 
Gaussian probability distribution and corresponds to 
taking only the N° terms in Eq. m, that is, only 
the k = 1 terms in the first two sums on the right- 
hand side of this equation. This is displayed using a 
dashed line (purple online). The dot-dashed (blue) 
line and the solid (black) line refer respectively to 
keeping terms up to and including 1/N 3 ^ 2 and l/N 2 , 
respectively, in Eq. m The (green) diamonds rep¬ 
resent the distribution rebuilt from direct stochastic 
simulations, based on the Gillespie’s algorithm, for 
the system. When the order of the approximation is 
increased, the theoretical distribution tends to adjust 
well to the numerical profile. 

In Fig. [2] the performance of the van Kampen 
system-size expansion and the WKB scheme are com¬ 
pared with each other and with the results of simula¬ 
tions. The agreement between the two schemes and 
the numerics is satisfactory, with the overall skewness 
of the distribution appearing to be correctly captured 
by the approximations. The WKB solution tends to 
deviate from the expected profile for negative values 
of the fluctuation £ (at system-sizes of a few tens of 
individuals), while the van Kampen approximation, 
at the order of the expansion that we are working at, 
still proves to be adequate. The fact that the van 
Kampen scheme matches the simulated data, consti¬ 
tutes an a posteriori validation of the closure strategy 
implemented. 



Figure 1: The (green) diamonds represent the distribution 
in lin-log scale rebuilt from Gillespie’s algorithm simula¬ 
tion. The (purple) dashed line, the (blue) dot-dashed 
line and the (black) solid line show the predicted profile 
obtained with the first thirty moments and Eq. 113 trun¬ 
cated at order one, 1/A 3 ^ 2 and l/N 2 , respectively. The 
parameters used are: N = 1000 c = 0.5 b = 0.3 d = 
0.2. For this choice of the parameters the fixed point of 
the deterministic dynamics is <j >* = 0.125. 



Figure 2: The (green) diamonds represent the distribution 
in lin-log scale rebuilt from Gillespie’s algorithm simula¬ 
tion. The (black) solid line represents the van Kampen 
approximation with the first thirty moments, and the ex¬ 
pansion taken to order of l/N 2 . The (red) dashed line 
is the WKB approximation. The parameters used are: 
N = 1000 c = 0.5 b = 0.3 d = 0.2. For this choice 
of the parameters the fixed point of the deterministic dy¬ 
namics is (j>* = 0.125. 


IV. CONCLUSION 

In this paper we have investigated the distribution 
of fluctuations in a stochastic model which was cho¬ 
sen to give the logistic equation as its deterministic 
limit. This equation has the necessary features of a 
stable non-trivial fixed point and an unstable trivial 
fixed point required to explore extinction dynamics. 
Stochastic fluctuations then take the system from the 
vicinity of the non-trivial fixed point of the determin- 
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istic dynamics, to the trivial fixed point, where all 
the individuals have become extinct. For sufficiently 
large population sizes, extinctions will be rare and 
the state of the system will fluctuate about the non¬ 
trivial fixed point for long time. 

In this situation the van Kampen system-size ex¬ 
pansion constitutes a powerful analytical tool to es¬ 
timate the distribution of fluctuations. Essentially 
the method gives the deterministic equations to lead¬ 
ing order in an expansion in (inverse) system-size. 
The linear stochastic corrections found at the next- 
to-leading order correspond to Gaussian fluctuations, 
an excellent approximation to the dynamics if ex¬ 
tinctions are so rare as to be negligible. When the 
size of the population is reduced, non-Gaussian traits 
prove crucial and one needs to go beyond this conven¬ 
tional order of approximation to eventually resolve 
the skewness of the distribution. 

Alternatively, a WKB-like perturbation scheme 
can be implemented to derive a closed analytical ap¬ 
proximation for the distribution of fluctuations. The 
method consists of postulating that the dominant 
contribution to the probability distribution behaves 
for large system size, A, as e ~ NS ( x ) where S(x) is 
a solution of a Hamilton-Jacobi equation found from 
the carrying out the WKB analysis. The correspond¬ 
ing Hamiltonian can be used to describe the extinc¬ 
tion trajectories, and consequently the distribution of 
stochastic fluctuations. 

The aim of this paper is to explore the connec¬ 
tion between the van Kampen system-size expansion 
at higher order and the WKB method. The two 
methods have a different basis, and it is therefore in¬ 
teresting to assess their respective predictive ability 
in regimes where extinctions are important. To this 
end, we have carried out the explicit calculations in¬ 
volved in these two methods, for the specific case of a 
stochastic model of the logistic type described above. 
The theoretical predictions have been compared to 
the results of the numerical simulations, obtained by 
solving the relevant stochastic model via the standard 
Gillespie algorithm 0- In both cases the agreement 
is satisfactory. 

The advantage of the WKB method over the gen¬ 
eralized van Kampen expansion is that the former 
enables one to recover a closed analytical expression 
for the distribution of fluctuations. In contrast, the 
latter requires working with a large set of algebraic 
equations for the moments of the distribution, a step 
that can be only performed numerically. In addition, 
the extended van Kampen method must be accom¬ 
panied by a dedicated truncation strategy, to get a 
fully consistent set of equations for the unknown mo¬ 
ments. The validity of this closure can only be tested 
a posteriori by a direct comparison with numerical, 
or experimentally available, data. On the other hand 
the WKB approximation is unable to capture the be¬ 
haviour near the extinction boundary; to do this it 


would need to be matched to a boundary-layer solu¬ 
tion calculated outside of the WKB scheme [19, 20 . 
However there is no need to do this in the cur¬ 
rent case, since the van Kampen expansion taken to 
higher-order fulfills this role. Thus, when used in con¬ 
junction with one another, the two approaches show 
that the extinction dynamics can be correctly cap¬ 
tured. 


Appendix A: The system-size expansion to 
higher orders 


In this Appendix we give some of the intermediate 
steps in the derivation of the van Kampen system-size 
expansion to all orders in A -1 / 2 . 

The starting point is the substitution of the Tay¬ 
lor series expansion given by Eq. © into the mas¬ 
ter equation (fol) . while using the van Kampen ansatz 
Eq. ©. For example, the first term in the master 
equation gives 

OO h -| 

(£+ - 1) [T(n - l|n)P„«)] = (Ep^»l) 

k= 1 ' ^ 

{d(j) + d^/VN+ 0^+ 2c^/VN+ ci 2 /N)H(i) . 


It is convenient to group these terms according to the 
power of £ which multiplies n, as follows: 


k=1 S 

OO 1 1 

+ ( fc! JSfk/2+1/2 Q£k ( 2c ^ + rf £) n ) 
k— 1 ' S 


/ ” i i d k 

V ■' k\ _/V fc / 2 + 1 

k= 1 s 



(Al) 


The lowest order contribution is the k = 1 term 
in the first sum, and this matches the second term 
in Eq. J7J) (after the time-rescaling to introduce r). 
Omitting this term in the first sum therefore gives the 
contribution to dH/dr from (£ + — l)[T(n—l\n)P n (t)\. 
Carrying out the same steps for the other term in 
Eq. ©, and taking out a factor of A -1 from the 
rescaling, gives the following expression for dH/dr: 


o° 

J2^Nkj2^T [^ + ^ 2 + (-l) fc (&^-&</> 2 )] 


k=2 


oo 


1 


k=1 


k\ Nk/2 -1/2 


[2c</> + d + (—l) fc (6 — 2b(f>)] 


fc=l 


d£ k 


d k U 

+ 

d k [£H] 

d£k 

(A2) 
















If we now introduce the simplifying notation 


we obtain the generalized Fokker-Planck equation 
HD of the main text. 


fk{4>) = [# + C( t> 2 + - be/) 2 )] , 

gk(4) = [2c</> + d + (—l) k (b — 2b(f>)] , 

h k = [c-b{- l) fc ], (A3) 
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